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t^ i Abstract. A fast harmonic oscillator is linearly coupled with a system of Ising spins 

that are in contact with a thermal bath, and evolve under a slow Glauber dynamics 

at dimensionless temperature 9. The spins have a coupling constant proportional to 

the oscillator position. The oscillator-spin interaction produces a second order phase 

' ^ '_ transition at 6 = I with the oscillator position as its order parameter: the equilibrium 

position is zero for > 1 and non-zero for 9 < I. For 9 < 1, the dynamics of this 
^ , system is quite different from relaxation to equilibrium. For most initial conditions, the 

oscillator position performs modulated oscillations about one of the stable equilibrium 

positions with a long relaxation time. For random initial conditions and a sufhciently 

^ I large spin system, the unstable zero position of the oscillator is stabilized after a 

^S^ ' relaxation time proportional to 9. If the spin system is smaller, the situation is the 

QQ , same until the oscillator position is close to zero, then it crosses over to a neighborhood 

of a stable equilibrium position about which keeps oscillating for an exponentially long 
f'*^. ■ relaxation time. These results of stochastic simulations are predicted by modulation 

^3 I equations obtained from a multiple scale analysis of macroscopic equations. 
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1. Introduction 

Many physical processes are modeled by means of an oscillator coupled to a thermal bath 
or to spin systems. Among numerous examples, we can mention the classical version of 
the spin-phonon system describing the collective Jahn- Teller effect [U |2] , decoherence of 
a spin representing a two-level system due to coupling to a boson bath (the spin-boson 
system) [3], a classical oscillator coupled to a spin causes wave function collapse thereof 
[1] , single molecule magnets or nuclear spins modeled as large spins coupled to a boson 
bath [5|, mass spectrometry through a nanoelectromechanical oscillator whose resonant 
frequency decreases as single molecules are added thereto [6], a 1/2-spin representing 
a nonlinear Josephson phase quantum bit is coupled to an oscillator (superconducting 
resonator) and to a classical signal [HIH], etc. 

Recently, we have studied a mechanical system (a classical harmonic oscillator) 
coupled to a chain of Ising spins in contact with a thermal bath at temperature T, so 
that 

* 1=1 

The spins have an energy —ixqX ^j^j^ a^aij^i/ \/N , and they flip stochastically according 
to Glauber dynamics at a rate Wi{cr\x,p) = (a/2) [1 — 7(Ji(crj_i -|- crj+i)/2], with 
7 = tanh.{2 fXQx/{kBT\/N)) [9], [10]. There is a second order phase transition at a critical 
temperature Tc = jJ^limojQkB) and the equilibrium position of the oscillator is its order 
parameter. Above the critical temperature, the oscillator equilibrium position is the 
same as that of the uncoupled oscillator. Below the critical temperature, two symmetric 
nonzero equilibrium positions issue forth from zero as in the diagram of a pitchfork 
bifurcation. In the limit of fast relaxation of the spins compared to the natural period 
of the oscillator, uo/a <^ 1, and ignoring fluctuations, the oscillator position satisfies 
an effective equation having a nonlinear friction and a nonlinear force term [10]. The 
stationary solutions of the reduced dynamics coincide with the oscillator equilibrium 
positions for each temperature. Numerical simulations confirm the theory except for 
very low temperatures (sufficiently lower than T^). In this temperature range, the 
numerical simulations show underdamped oscillations towards a nonzero equilibrium 
position of the oscillator while the theory predicts monotonic evolution towards the 
same value pLQj. Breakdown of the theoretical predictions is expected for sufficiently 
low temperatures because the relaxation time of the spins increases as the temperature 
decreases [HI [121 [13 [H] ^^^ the assumption that the spins relax in a time scale much 
shorter than the oscillator natural period is no longer true. 

In this paper, we consider the same model of an oscillator coupled to Ising spins 
that evolve according to Glauber dynamics but in the opposite limit of a short oscillator 
period compared to the spin relaxation time, i.e. ujq <^ a. The main characteristics 
of the model are reviewed in section [2l We use nondimensional variables in which the 
natural frequency of the oscillator and the transition temperature are both equal to one. 
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In a mesoscopic description, the dynamical equations of the system are stochastic due to 
the couphng of the spins with the thermal bath at temperature T. These equations are 
equivalent to a master equation for the joint probability density V{x,p, (T,t) of finding 
the oscillator with position x and momentum p, and the spins in configuration cr at time 
t. Its stationary solution is the canonical distribution, which describes the equilibrium 
behavior of the system, including the second order phase transition described in the 
previous paragraph. The proof showing that it is the only stationary solution and that 
any solution of the master equation for finitely many spins tends to it in the long time 
limit is outlined in Appendix A 



The rest of the paper is as follows. In section [2], we recall the model introduced in 
[To] and derive the nondimensional macroscopic equations valid in the limit of infinitely 
many spins and in the absence of fluctuations. In section [HI we analyze these equations 
in the limit of large oscillator frequency as compared with the reciprocal relaxation time 
of the spin system, S = a/uo <C 1. We show that there are two well separated time scales 
governing the behavior of our system: the fast time scale corresponding to the natural 
period of the oscillator and the slow time scale t = 6t over which the oscillator and 
the spin correlations relax to equilibrium. The oscillator position x^^^ is approximately 
given by the order one spin correlation, Ci = A^""*^ X]j=i('^«^i+i) (^^ ^ ~^ ^^)> P^^^ ^ 
modulated oscillatory term. The equations for the spin correlations are the typical ones 
for Glauber dynamics, but with their parameters averaged over an oscillator period. 

The main results in this paper are summarized in section [31 The averaged equations 
for the slowly varying spin correlations and the envelope of the fast oscillations in x^'^^ 
constitute our first result. These averaged equations are approximated in two limits 
corresponding to our results 2 and 3: near the critical temperature 6 = 1 and at 
low temperatures 6-^0. For most initial configurations producing a nonzero Ci, the 
correlations evolve towards their equilibrium values. Near the critical temperature there 
is critical slowing down (on a time scale |1 — 0|r) in the already slow evolution of Ci 
toward its equilibrium value while the oscillation envelope decays to zero over the time 
scale T. The behavior for low temperatures is more surprising. The correlations decay 
algebraically (long time tails) to equilibrium, while the amplitude of the oscillations 
in the oscillator position does not vanish (except on an exponentially long time scale 
which is outside the scope of standard numerical calculations). On the other hand, 
if the initial configuration of the spin system is random, the system approaches the 
unstable equilibrium solution Ci = 0, 5; = for ^ < 1 (we have checked this for 
temperatures as low as 6' = 0.37) instead of going to the stable equilibrium solution 
with nonvanishing Ci and x. For initial conditions arbitrarily close to the random ones, 
the system first approaches the unstable equihbrium solution until it is quite close to it. 
Then it crosses over towards the stable solution having nonzero Ci while the oscillator 
position undergoes fast oscillations about Ci with an exponentially small damping. This 
departure of the equilibrium behavior for exceedingly long time intervals is reminiscent 
of dynamical glassy behavior in spin glasses and other slow relaxing systems [T5l [16] . 

The comparison between our results and numerical simulations is excellent, as 
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discussed in section HI Section \5\ contains detailed derivations of the results given 
in section |3l Section [6] contains our conclusions and the appendices are devoted 



to proving the H theorem for our system (Appendix A), and to technical details 



on the approximations used near the critical temperature (Appendix B) and for low 



temperatures (Appendix C). 



2. Dynamics 

Let us summarize the main aspects of the model introduced in [10] , where more details 
can be found. The system consists of a one dimensional harmonic oscillator (mass m, 
frequency cjq, position x and momentum p) and A^ ^ 1 internal degrees of freedom 
modeled by Ising spins {ai = ±1, i = 1,...,N) in contact with a heat bath at 
temperature T. The system has an energy 

2 1 ^ 

nj( \ P , ^ 2 2 /^O V;^ /r, \ 

* j=i 

The first two terms on the right hand side (rhs) of f l2aj) correspond to the energy 
of the uncoupled oscillator, while the last term stands for the interaction energy 
between the oscillator and the spins. The latter can also be understood as a nearest 
neighbour interaction between the spins with a coupling constant Jes = f^ox/vN which 
is proportional to the oscillator position x. The parameter /ip measures the strength of 
the coupling between the oscillator and the Ising system. 

The dynamics of the oscillator is governed by Hamilton's equations of motion, i.e. 

N 

X + WqX = ^ V (Tia.+i , (3) 

my/N "j-^ 

and the Ising spins evolve according to Glauber spin flip dynamics at temperature T. 
Thus at any time t, the system may experience a transition from {x,p, cr) to {x,p, Ricr) 
at a rate given by [9] 



a 



7 



2/ioX 



Wi{(T\x,p) = - 1 - -(Ti(ai_i + a,+i) , 7 = tanh , (4) 

Z I Z J \kb-I. V A/ 

where RiCr is the configuration obtained from cr by flipping the z-th spin, ks is the 

Boltzmann constant and T is the temperature of the system. The parameter a gives 

the characteristic attempt rate for the transitions in the Ising system. 



X p V Wi t 



Table 1. Nondimensional units and parameters. 

It is convenient to introduce nondimensional variables according to x* = x/[x], 
- t/[t], ..., where the units [x], [t], etc are as defined in Table [TJ Dropping the 
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asterisks so as not to clutter our formulas, we obtain the nondimensional equations: 



1^ 



+ X 



1 ^ 



(5) 



1 A 



^t + P^^+ \jt2_^ o-iO-j+i 



X dr, 



i=l 



V{x,p,(T,t) 



N 



6Y,[Wi{RifT\x,p)V{x,p,Ri(T,t)-Wi{(T\x,p)V{x,p,(T,t)], (6) 



i=l 



Wi{cr\x,p) 



7(x) 



0-i[(Ti_i +0-J+1 



7(x) 



, ,'2x 
tann ( — - 
9 



Here 6 and 6* are dimensionless parameters given by 
a -^ T_ y^_ /ig 



5 



9 



(7) 
(8) 

(9) 

5 is the ratio of the characteristic spin rate a to the oscillator natural frequency ojq 
and 9 is the dimensionless temperature. At ^ = 1 (T = Tc), there is a second order 
phase transition whose order parameter is the equilibrium position of the oscillator 
[in]. Equation (jS]) is equivalent to the oscillator equation ([3]), while (jH]) is the master 
equation for the joint probability density V{x,p, cr, t) of finding at time t the oscillator 
with nondimensional position x and momentum p, and the spins in a configuration 
(T = {cTi, cr2, . . . , cttv}. The equilibrium canonical density is the stationary solution 



Veq{x,p,(T) 



n{x,p,cr) 



4^-p 


-j'H{x,p,(t) 


p' + x' 


N 


2 



(10) 

(11) 



i=l 



where the partition function Z guarantees that Vcq{x,p,cr) is normalized to unity. 
Summing over all the spin configurations cr, we derive the marginal probability Vcq{x,p) 



^eq(a;,p) = -exp 



A^ fp 
T 



+ V,s(x] 



Veflf(x) 



2 



In cosh ( — ) + In 2 



:i2) 



(13) 



For finite N, the H-theorem of [Appendix A| proves that ( jTOl) is globally stable and that 
it is the only stationary solution of the master equation iQ. The maxima of Veq{x,p), 
(xeq,Peq), coiucidc with the equilibrium mean values of x and p in the limit as A^ — )■ oo. 
They are given by the extrema of Veff(a;) in (TT3|) . i.e. by the solutions of the equation 



dx 



0. 



X. 



eq 



ceq 



0. 
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Clearly Xgq = is always a solution for any 9. For ^ > 1, it is the only solution, it 
corresponds to a maximum of Pcq and is therefore stable. Ai 9 = 1 two new stable 
equilibria bifurcate from Xgq = and exist for 9 <1. As^— >l^,we have 

Xeq~±V3(l-e). (15) 

In the limit as A^ — > oo, we split the variables x = x + Ax, 0"jCrj+„ = C„ + ACj^„, 
where x = (x) and C„ = {aiai+n) (independent of i, provided we only consider 
homogeneous initial correlations {cri(yi+n){t = 0) that depend on n but not on i). Then 
we insert the result in the following equations obtained from ([S])-(IHD, 

^ + w = ^f;(a,>, (16) 

i=l 
1 r 

^^{C^,n) = ^{l{x) {Ctn-l + Q,n+1 + Ci-l,n+l + Q+i^n-l)) - 26{a^n) (17) 

(for n > 1 and 2 = 1,..., A^), and ignore the fluctuations Ax and ACj^„. We obtain the 
following system of equations for the macroscopic variables x and C„ [lOj : 

d^T — 

^ + ^ = Cu (18) 

j^Cn = 6^{x)(c^i + a^i^-26Cn,n>l, Co = 1, (19) 

7(x) = tanh (^] . (20) 

The mean-field or macroscopic dynamical behavior of the oscillator-spin system is found 
by solving the equations (IT5]1 and (IT^ with the boundary condition Cq = 1 and 
appropriate initial conditions. 

3. Fast oscillator dynamics. General results 

The limit of a fast oscillator compared to spin relaxation describes the region of low 
temperatures of our system (no matter what the oscillator natural frequency is) because 
the spin relaxation time becomes arbitrarily large as the temperature decreases towards 
zero. We have obtained the following results: 

Result 1. In the limit as 5 — )■ 0, the system of nondimensional equations ( fi^) - ( fTPj) 
can he approximated by the multiscale solution: 

xit;6)=x^'\x,r) + Oi6), a;(t; <5) = Cf (r) + 0(5), (21a) 

x(°)(x,r) = cf(r)+i?(r)sinx, r = <5t, x = t + 0(r), (216) 

in which C„ (r), R{t) and 0(r) solve 

^ = lixi0Kx,r))iCl^\ + d:i)-2Ci^\ Cr = l, (22a) 

^ = -(l + Cf)7(xW(x,r))sinx, ^ = 0. (226) 



Nonequilibrium dynamics of a fast oscillator coupled to Glauber spins 7 

Here we have defined the time averages over the fast periodic variable x (keeping the 
slow variable r fixed) as 

7b() = ^ r f{x)dx. (23) 

This result will be derived later in section |H Note that, according to (121 gp and 
(12 1 b\i . the oscillator velocity is 

^^x{t;6) = ^{t + <p,6t) + 0{6). (24) 

Then the initial conditions 

xiO;5) = Xo, i{0;5) = vo, 5^(0; 5) = ^^(O) (25) 

yield the following initial conditions for the approximate quantities: 



R{0) = ^ [xo - C,{OW + vl sin0=^^-^^, (26) 

for R{0) i- 0. If Xo = C^(0) and Wq = 0, i?(0) = and ^^ produces R(t) = so that 
there are no modulated oscillations for this choice of initial data. = f o = dx/dt = 
yields dCf^/dr = at t = 0. 

The stationary solutions of (122 ap and (122 6p . Cn,eq and Xcq , are constants, 



independent on x ^"^^ t- Therefore i? = 0, xiq = C*i,eq5 ^^^ '^^ have 



m , r^(o) 



2^n,eq + 7(^l,oq) ^n-l,eq + ^n+l,eq ~ U , ^0,eq ~ ^ • \^') 



m 



This system has the unique solution 



with 



Cn,eq " C';^ ) , (28) 



m 



which determines the equilibrium value of C| or x^*^-*, since both of them are equal. 
Equation fl2^ is equivalent to dH]), the bifurcation equation. Then the stationary 
solutions of the approximate dynamical equations (I22ap - (]226p coincide with the exact 
equilibrium solutions obtained in the previous section. There is only one stationary 
solution R = Cn = {n > 0) ior 6 > 1 (which is stable) and three stationary solutions 
i? = 0, C^iq = (cJJqV, with C[% solving ([29]) for < 1. In the latter case, CJ^^ = 

is unstable whereas the solutions with nonzero CJ j are stable and bifurcate from the 
zero solution at 6' = 1. This was to be expected, since the equilibrium distribution does 
not depend on 6, and it shows the consistency of the multiple scales approximation. 

Equations (\22a!} - (\22b\} have solutions Cn = {n > 1), with R{t) decreasing from 
i?(0) > to as r — 7- oo {dR/dr < because 7(i?(r) sinx) sinx > in x G [— 7r,7r]. 
Then its average appearing in (\22b\} is positive for i? > and it vanishes for R = 0). 
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8 



For ^ < 1, this solution lies on the stable invariant manifold of the saddle point R = 0, 
Cn = (n > 1): linearizing fl22ap - fl22&|) for n = 1 about the saddle point, we find 



dC| 



(0) 



1 



2(^-l)C^l 



(0) 



(30 a) 
(306) 



dr 
dR _i? 

d7 ^ ~J' 

Setting C„ = (n > 1) in (IT8|) -( |20|) we no longer obtain an invariant manifold as it is 
the case with Cn = for the averaged equations {\2'2a^ - (\22b^ . In fact, ( 1T^ becomes 
dCi/dt = 7(x) ^ for n = 1 if C^ = (n > 1) and O) = 1. 

There are two limits in which we can find approximate solutions of the averaged 
equations (122 ap and ( 122 6^ . namely 6* — )• 1 (the critical temperature) and 6* — )■ 0. 



3.1. The limit e ^ 1 



Result 2. In the limit as 6 -^ 1, the solutions Cn (t) and R{t) of \22a\) and Ii22b\) are 
given by 



m 



Crir) = [signCrm 



m, 



\ 



3(1 



1 _L p-4(l-e)r ( 3(i-e) 



0(|1-^|), (31a) 



(316) 
(31c) 



Ci'\T) = [C?\r)ni + 0{\l-e\'/')],n>2 
R{T) = R{0)e-^ + O{\l-e\), 

provided the initial conditions -R(O) and C| (0) are both of order e = a/|1 — d\- 

These formulas will be derived later in section H] by using bifurcation theory (we 
understand sign(O) = in the previous formula, so that CJ (r) = if C| (0) = 0). In 
Appendix B , we show that ( 131 aP and ( 131 6p with r = 6t also hold for Ci = x and C„, 
respectively, except for exponentially decreasing terms. 



3.2. The limit e -^0 



Result 3. In the limit as 6 ^ 0^ , Equations ^22a\) and ^22b\) can be approximated by 



dCf 

dr 
dR _ 

d7~ 



- D \A-yn-l + ^n+lj ~ '^^n i ^0 

-2(l + Cr)e-'l^ri/»A(^). 



.{0)| 



.(0): 



with S =signC{ (t), whenever R < \C\ \, and for R > \C\ \ by 



dr 



l{d:i, + Ci^i:)^rcsin(^] 






dr TT 



\ 



CI 



(0)' 



R 



(0) 



(32a) 
(326) 

(33a) 
(336) 
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Note that rescaling Cn -^ S''^Cn , we can reduce the case S' = — 1 in ( 132 aP to the 
case S = 1. For C| < 0, i? and Cn for even n remain the same as in the case C| > 
whereas the correlations Cn for odd n change their sign. For CJ > 0, 

oo 

Cf (r) = 1 - e-'- 5^[1 - Cf(0)] [/,„_, I (2r) - /„^.,(2r)] , (34) 

i=i 
where In{z) is the modified Bessel function of the first kind with index n, solves (132 ap . 
Then the asymptotic properties of In{z) [H] yield as r — > oo the long time tail result 

oo 

C(°)(r) ~ 1 - n(4vrr^)-V2^,[l - cf (0)], (35) 

provided that the sum is finite. As r — > oo, i \32b\i becomes 

f^-4e-.».(f), (3a, 

i? — )■ and, using the relation /i(x) ~ x/2 [T7] in (15^ . this equation becomes 

^^_^e-^/^i? (r^oo). (37) 

R decays exponentially with a characteristic time which diverges exponentially as ^ — t- 0. 
The amplitude of the oscillations around the steady value x^q = 1 is damped only after 
an exponentially large time ^e^/^/8. 

Let us consider now equilibrium-like initial conditions for the Ising system, C„(0) = 
r", with < r < 1. In this case, the solution flMl) becomes 



Cl^\r) = ^ fdg ;f f e-^-(^— ^) + Te"- [/„_,(2t) - 4^,(2t)] dt. (38) 
TT Jo 1 + r^ - 2r cos q Jq 

For n = 1, fl38|) increases monotonically from CJ (0) = r to its asymptotic value 1. 
Since R decreases with time, C| (r) > R{t) if this condition holds at r = 0. Using the 
asymptotic expression of In±i{2t) for large t [IT], we find 

Cf)(r)~l--^ + 0(r-3/2). (39) 



'vrr 

For this subset of initial conditions, the long time behavior of the Ci is independent of 
the initial conditions. It must be noted that this is consistent with (1351) . since for these 
initial conditions the sum in fl35|) diverges. 

Let us now consider equilibrium-like initial conditions with < r <^ 1, so that 
i?(0) > C[°\0) > 0. Then ([SSaD-dHMI become 

2^ = ^C™-2Cr, Cr = l, (40a) 

^ ^ -'- (40.) 

clr n 

We have assumed that C^°^ = 0(r") and i?(0) = 0(1) as r ^ 0. Then (gOoD yields 
^(0) ^ ^(o)^Q^|^_(l_2/^)r g^^^ gog gj^gg j^(^^^ ^ j^(^Q^ _ 2^/^^ 'ji^g correlations remain 

small whereas R{t) becomes zero at a time 



ro = ^ = l\/[xo-CM? + vl (41) 
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At a slightly smaller time than this, R{t) ~ C\ (r) and the approximations ( 140 ap 
and ( 14061) break down after some time proportional to ln(l/r). If r = 0, then we 
are on the stable invariant manifold of the stationary saddle point solution mentioned 
above, all correlations remain zero and R will decay to zero at r = tq and remain there 
indefinitely. However the stationary state R{t) = CJ (r) = is unstable as indicated 
by the linearized equation ( 130 al) and any numerical error could give rise to a "virtual" 
small r and lead to the neighborhood of one of the stable stationary (equilibrium) 
solutions after the break down time before mentioned. The solutions of the macroscopic 
system (IT8l) -( l20l) will exhibit a different behavior because this system does not have the 
same invariant manifold as the averaged system (I22al) - (I226|) . These conclusions will be 
checked with numerical solutions in section HI 

4. Comparison with numerical simulations 

To check the validity of our approximations, we shall first compare the solutions of the 
macroscopic equations (IT8l) - (IT9l) to the averaged system (I22ap - (I22?)|) (Result 1 in section 
E]) and to one of the approximations thereof: either bifurcation theory for 9 close to 1 
(Result 2, equations ( l31ap -( l3Tcj) in section [3]) or the set of two equations (I32ap -( l33^ for 



low temperature 9 (Result 3 in section [3]). Later we will compare our theory with direct 
numerical simulations of the oscillator-spin system. 

4-1. Comparison with numerical solutions of the macroscopic equations 

Figures [H and [2] compare the direct numerical solutions of the macroscopic equations 
to the approximations given in Results 1, 2 and 3 of section |3l Near the critical 
temperature, figure [1] shows that the approximations given by Results 1 and 2 (averaged 
equations and bifurcation approximation) are excellent for sufficiently small 5 (0.1 
and smaller values), whereas there is a small transient that is not captured by the 
approximations for larger values of 5. For those larger values, the approximation 
of the averaged equations given by bifurcation theory (Result 2) is somewhat better 
than the averaged equations themselves (Result 1). This is not surprising as the 
bifurcation equation for the averaged equations (I22ap - (l226p is the same as that for 



the full macroscopic equations as shown in [Appendix B[ An even better approximation 



can be found if the initial time layer is built by directly using the transients for the 
macroscopic equations. For low temperature, 9 = 0.37 (|xoq| = jC'i^eql — 0.99), figure 
|2] compares the solution of the macroscopic equations to the approximations given by 
Results 1 and 3 in section |3l We observe that even the rougher approximation given by 
Result 3 ranks from good to excellent as 6 decreases. 

Next, we probe the dependence of our approximation on initial conditions. In figure 
ini we fix = 0.1 and 6 = 0.01, and use equilibrium-like initial correlations Ci(0) = r" 
for different values of r that range from very small to near 1. We find that Results 1 
and 3 approximate very well the solution of the macroscopic equations except for very 
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Figure 1. Numerical solution of the macroscopic system (solid red line), Result 1 
(dashed blue line) and Result 2 (dot-dashed green line) for the following values of 5: 
(a) 1, (b) 0.5, (c) 0.2, (d) 0.1, (e) 0.05,^) 0.02, (g) 0.01. J^n aU cases, 9 = 0.98, 
and initial data are xq = 0.1, wq = 0.1, Ci(0) = xq + 0.1, C„(0) = 0, n > 2. The 
macroscopic system has been truncated by imposing Ciooo = 0- 
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Figure 2. Numerical solution of the macroscopic system (solid red line), Result 1 
(dashed blue line) and Result 3 (dot-dashed black line) for the following values of S: 
(a) 1, (b) 0.5, (c) 0.2, (d)^0.1, (e) 0.05, (f)J).02, (g) 0.01. Here = 0.37, initial data 
are xo = 0.5, vq — 0.1, Ci(0) ~ xq + 0.5, C„(0) = 0, n > 2, and the other parameter 
values as in figure [TJ 
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small r. As r decreases, the system first tends towards the unstable solution x = 
(attempting to stabilize it) and then it crosses over to the basin of attraction of the 
stable equilibrium with x close to 1. For r = (or r ^ 1), the approximations given by 
Results 1 and 3 indicate stabilization of x = at the time ( HTj) . i.e. to = '''o/^ = 80.09. 
This is confirmed by figure M,g) and it agrees with the fact that Cn = (n > 1), 
R = R{t) is an invariant stable manifold of the averaged equations (I22ap - fl22?)|) . For the 
case r = 0, R and Cn should remain zero. However, figures [3]^g) and (h) show that 
numerical errors build up and send R{t) given by the solution of the averaged equations 
(\22a\i - l\22b\} to one of the stable equilibrium values (in the case of figures [3]^g) and (h) 
to the wrong one, near x = —1) after some time ti > tq. Since zero correlations do not 
correspond to an invariant stable manifold of the macroscopic equations flT8|) - fl20|) . the 
solutions thereof depart in the vicinity of zero and after some time from the solutions 
of the averaged system (I22al) -( l226l) . cf figures [3](g) -(h). 

4-2. Comparison with numerical simulations of the oscillator-spin system 

We now compare our approximate theories to direct numerical simulations of the 
oscillator-spin system with Glauber dynamics. We choose a; = 10 and a = 1, so that 
6 = 0.1. We have carried out simulations for a temperature 6 = 0.995, very close to 
the critical temperature (figure 111(a)), which corresponds to e = VO.OOS = 0.07. A large 
number of particles must be considered in the simulations, because of the divergence of 
the fluctuations at the critical temperature. Namely, the number of spins is A^ = 10^ and 
we have averaged over N^ = 1000 trajectories. The initial conditions have been chosen 
such that xq = Ci(0) and fo = 0. Therefore, we expect that -R(t) = for all times, and 
x{t) = Ci{t), both of them relaxing to their equilibrium value Xeq = C*i,eq — 0.122. Two 
different initial values of x(0) have been considered, one above the equilibrium value 
Xq = 0.2, and one below it, xq = 0.1. In both cases, the simulation curves agree with 
the theoretical prediction given by (131 ap . For a lower temperature, 6 = 0.98 the theory 
still gives a good description of the simulation results while the initial conditions are 
not too big (figure m^b)). This value of 9 corresponds to e = VO.02 = 0.14, which is not 
so small (in fact, the steady value for the oscillator position is 0.24, and its maximum 
value is unity). A further test of the theory is done in figure [5l for the same values of 
the parameters as in figure Hl^a), but xq 7^ Ci(0). The initial conditions are Ci(0) = 0.2 
but Xo = 0.1 with vanishing initial velocity, so that we obtain from ( l26l) 



R{t) = O.le"^ , = -- ^ x{t) - Ci{t) = -O.le"^ cost . (42) 

In figure [5], we observe that the simulations of the macroscopic equations, the averaged 
equations (Result 1) and bifurcation theory (Result 2) all follow exactly (H2]) and agree 
very well with direct stochastic simulations. 

Considering now equilibrium-like initial conditions for the Ising system, i.e. C„(0) = 
r" with < r < 1, the low temperature correlations are given by ( 138|) . In figure Ej we 
have chosen parameter values 6 = 0.37, S = 0.1 and initial conditions Xq = 0.25, fo = 0, 
r = 0.2. In this way, as xq — Ci(0) is nonzero, we expect to observe the damped 
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Figure 3. Numerical solution of the macroscopic system (solid red line), Result 1 
(dashed blue line) and Result 3 (dot-dashed black line) for 9 = 0.1, S = 0.01, xq = 0.5, 
vo = 0.1, Ci{0) = r" with r =: (a) 0.9, (b) 0.5, (c) 0.1, (d) 0.01, (e) 0.001, (f) 
same as (e) but including results of stochastic simulations (magenta circles), (g) 0. (h) 
corresponds to S — 0.05 and r — 0. Note that the approximation given by (|33aP - ()33&P 
(Result 3) breaks down when R — and the corresponding line is interrupted in (g)-(h) 
after the breakdown time. 
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Figure 4. Averaged trajectories {x(t)) = x{t) for: (a) = 0.995, vq = 0, Nt = 10^, 
Ci(0) = xq = 0.2 (upper curve, cyan squares) and Ci(0) = xq = 0.1 (lower curve, 
magenta circles); (b) 9 = 0.98, Nt = 10'^ and from top to bottom, Ci(0) — xq = 0.5 
(green triangles), 0.3 (cyan squares), 0.1 (magenta circles). Other parameter values 
are ^ = 0.1, iV = 10^. Solid red line: macroscopic equations ([T8|) - ([T9l) . dashed blue 
line: averaged equations (|22a|) - (j226p . dotted black line: bifurcation theory pia|) - pic|) . 
Most of these theoretical predictions are indistinguishable in the scale of the figure. 
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Figure 5. Difference {x(t)) - {Ci{t)) = x{t) - Ci{t) for Ci(0) = 0.2, xq = 0.1, vq = 0, 
9 = 0.995, S = 0.1, N = 10^, Nt — 100. Stochastic equations (magenta circles and 
cyan squares), macroscopic system (solid red line). Result 1 (dashed blue line) and 
Result 2 (dot-dashed black line). Also plotted x{t) - Ci(t) = for a;o = 0.2. All 
predictions are indistinguishable in the scale of the figure. 



oscillations of x{t) around Ci{t). The theoretical curve for Ci{t) is not plotted, because 
it is indistinguishable from the numerical one, and the oscillations of x{t) around Ci{t) 
are clearly observed. In figure [Tl^a), the difference x{t) — Ci{t) is shown, for the same 
values of the parameters as in figure |6l We observe clear modulated oscillations with an 
envelope R that decreases over a very long time scale, as shown in figure [Tl^b) (note that 
the time scale thereof is much longer than that in figure [Tl^a)). The long time behavior 
of R is very well described by the exponential decay in equation ( l37l) . 
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Figure 6. Averaged trajectories {x{t)) — x{t) (magenta) and (Ci(i)) 
for r = 0.2, xo ^ 0.25, vq^O^O^ 0.37, 5 = 0.1, N = 10^ Nt = 10^. 
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Figure 7. (a) Difference {x{t)) — {Ci{t)) = x{t) — Ci{t) from stoclrastic simulations 
with r = 0.2, xo = 0.25, wq = 0, 6* = 0.37, S = 0.1, TV = 10^ Nt = 10^. (b) Amplitude 
R (log-scale) vs time t. The dashed line is the best long time {t > 5000) exponential 
decayfit,i? = 0.045 exp(-4.39xl0"'5^). This is consistent with i? ex exp(-4.36xl0"6i) 
obtained from (|37p for very long times {t = 6t > 500). 



Next we consider the case of random initial conditions, C„(0) = for n > 1. 



-(0), 



For 9 < 1, fl22ap -( l226j) have the solutions C„ ' (r) = 0, n > 1, with R{t) decreasing 
monotonically to zero, according to (I30ap - (]306p . We have investigated this stabilization 
of the unstable equilibrium solution x^q = Ci,cq = by stochastic simulations of the 
oscillator-spin system, initially prepared in a completely random configuration. In 
figures [8] and [H we have chosen xq = 0.5 and Vq = 0.1, respectively, but we have 
observed the same behavior for all the other values of Xq and Vq we have tested. For 
6 = 0.95, figure [Ht^a) shows modulated oscillations of the oscillator position with an 
envelope that decreases monotonically to zero and remains there, in fact stabilizing 
the unstable solution Xeq = 0. We depict the oscillation envelope in[H](b), together with 
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the theoretically predicted exponential decay given by (130 5|) . The agreement is excellent 
(note the logarithmic scale in the vertical axis). It must be stressed that, in both figures, 
we are plotting only one trajectory for a very large system, A^ = 10^. All the individual 
trajectories (for fixed Xq and Vq) coincide. For a smaller system, the fiuctuations are 
larger and the system finishes in one of the two stable stationary states, Xeq > or 
Xeq < 0, depending on the trajectory considered. This is shown in M^c), where two 
different trajectories for A^ = 10^ are depicted. 

It is remarkable that the simulations of the oscillator-spin system for very large N 
resemble the separatrix trajectory of the averaged system (I22ap - (]226p . given by a R{t) 
that solves (12261) with CA = Sno- The averaged equations approximate the macroscopic 
equations (IT8|) - flT9|) obtained in the limit as A^ — ?■ oo by ignoring fiuctuations. However 
the separatrix of the averaged system that ends up at -R = is not a trajectory of 
( 1T8|) -(1T9l). so it is surprising that the approximate average equations describe better 
the result of stochastic simulations for random initial conditions than the more exact 
macroscopic equations. It also comes as a surprise that stochastic simulations for 
a smaller spin system with random initial conditions resemble the solutions of the 
macroscopic equations (TT8|) - (TT9l) for which the separatrix is not a trajectory. 




(a) 





Figure 8. (a) Oscillator position {x{t)) for r — 0, xo — 0.5, vq = 0, 9 = 0.95, 6 — 0.01, 
N — 10*, Nt = 1. This individual trajectory shows the stabilization of the unstable 
state Xeq = 0. (b) Amplitude R (log-scale) vs time, (c) Oscillator position {x{t)) for a 
smaller system, TV = 10^. We have depicted two typical trajectories, corresponding to 
different choices of the initial spin configuration. 
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We analyze the low temperature case {6 = 0.37) in figure O Parameter values 
are the same as in figure [H except for a larger initial oscillator position, xq = 2. This 
choice enlarges the initial time window during which R decreases linearly with time, cf 
( 140 b\i . Again, the oscillator approaches the unstable "trivial" solution and not one of 
the stable states with |a;eq| — 0.99; see figure |^a). The oscillation amplitude decays to 
zero: linearly with time according to ( \AOb\i in a first stage where 2R ^ 6, whereas it 
decays exponentially according to (130 6|) once 2R/9 ^ 1, see figure [9t^b). In figure [9]^c), 
we show the behavior for a smaller system with A^ = 10^. In this case, the fiuctuations 
are larger and they can drive the system out from the unstable state towards one of the 
two stable states, depending on the initial spin configuration. Again the simulations 
of the stochastic oscillator-spin system resembles the description given by the averaged 
equations for larger A^ and the description given by the "more exact" macroscopic 
equations (TT8l)-(TT9l) for a smaller A^. 




(a) 



(b) 




Figure 9. Same as in figure[8]for 6 = 0.37, xq — 2. In (b), R vs time shows the initial 
hnear regime, (J40&|1 . and the long time exponential decay, (|30&p . 



5. Derivation of the results in section [3] 

5.1. Derivation of Result 1 

To derive an approximate solution of (fT8|) - ( TT9|) as in Result 1 of section [3l we make 
the multiple scales Ansatz 

x{t; 6) = x(°)(t, r) + 6 x^^\t, r) + 0{6^), (43a) 
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C;(t;5)=Cf(t,r) + 5C«(t,r) + 0(52), (436) 

with r = 5t. While x and C„ are functions of one variable t, we consider that x'--'^ and 
Cn are functions of two independent variables, t and r = 5t. This is not really so, and 
therefore we are free to impose an additional condition to the multiscale functions in 
order to determine the evolution in the slow time scale r. This additional condition is 
that the two-time functions x'^^^ and Cn should be bounded for arbitrarily large values 
of t and any fixed value of r. We will see later that this condition indeed determines 
the slow-time evolution in the scale r. Insertion of ( 143 ap and (14361) in (ITSl) and (fT9l) 
produces the following hierarchy of equations 

^2^(0) 



dt^ 

dt 



-x^'^ = C['\ (44a) 

0, (446) 



/52 (1) f)2 (0) 

^^- ^{x^'W:\ + d:i)-2C^>^-^, (456) 



dt " '^ n-, ■ n+u n g^ 

and so on. 

Eq. (I2B1) implies that Cn is independent of t. Then the solution of (144 gp is 

a;(°)(x, r) = cf (r) + ^(^) sinx, X = ^ + 0(r)- (46) 

The oscillator position is therefore slowly driven by the spin correlation C| (r) while 
it performs modulated oscillations with slowly-varying amplitude R{t) and rapidly 
changing phase x- 

To find the slowly- varying quantities C„ (r), R{t) and 0(r), we analyze the first 
order equations (145 aP and (J456j). The rhs of (j456j) is 27r-periodic in the fast time t (or, 
equivalently, x)- Thus it can be expanded in a Fourier series in x whose zeroth harmonic 
should vanish for Cn to be bounded as t — t- oo. The zeroth order harmonic of the rhs 
of (H6l) is its average over a period of x- This implies Eq. (122 aP for Cn (t) with the 
definition (1231) for fast time averages. Eq. ( I22ap is similar to the Glauber dynamics with 
a corrected rate 7(x(°)(x, t)). Because of (H6|) . 7(x(°)(x, t)) depends on r through the 
functions R{t), (J){t) and CJ (r). 

To find -R(t) and 0(t), we should first integrate (145 6p . To this end, we write the 
periodic function 7(0;*^°^ (x, t)) as a Fourier series 

+00 
7(xW(x,r))= Y. 7.(r)e'^^ (47a) 

j=-oo 

7,(r) = ^ r dxe-'^'^7(a;^°nx,^)) = e-^x^(x(o)(x, r)) . (476) 



Substituting (122 aO in (145 6p . the latter equation can be written as 



dt 



7(a:(°))-7(a:(o)) C^l + ci"^ 



1+1 
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+00 
= E 7,(r)e^^->^(*-)[ci°2i(r)+Ci°i(r)] . (48) 

j=-oo,j^O 

This equation can be immediately integrated with the result 

00 , . 

CL'\t,T) = [d:\{T) + cZ{T)] Y. ^e^'^ + D^ir), (49a) 

j=-oo,j^O 



jj = e-'J>^7(x(°)) = 7(xW) cos(jx) - i 7(^(0)) sin(jx). (496) 

Here the "constant of integration" Dnir) is a slowly- varying function that can be 
determined from higher order equations in the hierarchy. Inserting (146|) and (149 ap in 
(I45ap . we find 






For x^^\x^ t) to be bounded as x — )■ 00, the right hand side of (150!) cannot contain first 
harmonic terms, and therefore the following equations hold: 

^ = -(1 + Cf)7(a:(o))sinx, (51a) 



.(0)^ 



i?-p = -(l + CD7(xW)cosx = 0, (516) 



which are (122 6p . Equivalently, Eq. (15 lap can be written as 



^ = ^(1 + Cf ) Ain (1 + e-4(cr+i^sin,)/.y (52) 

5.2. Derivation of Result 2 

At the critical temperature 6* = 1 we have a pitchfork bifurcation with nonzero stationary 
solutions that bifurcate from i? = C} = for 9 < I. We can obtain approximate 
solutions of the macroscopic system flTSj) and flT9|) directly by using bifurcation theory. 



See Appendix B However to keep the unity of our theoretical description, it is 
instructive to study the dynamical behavior of the reduced system (I22ap - fl226p near 
the critical temperature. The stationary state i? = C| = is stable for 9 > 1 and 
unstable for ^ < 1 so that the pitchfork bifurcation is supercritical: nonzero stable 
stationary branches exist for 9 < 1. From flTSjl . we anticipate that the stationary 
solutions bifurcating from x^^^ = Cn = are 

< = Cfl~±v/3ir^, 9^1-, (53) 

whereas the other correlations are given by fl28l) . Therefore, by defining 

9 = 1- e^92 , (54) 

a;eq and C*} j^ are of order e, while Ci,eq = 0(e"). The parameter e therefore measures 
the amplitude of the bifurcating solution. We have introduced the parameter 92 in ([5i 
to analyze simultaneously the cases < 1 (^2 = 1) and ^ > 1 (^2 = — 1)- 
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Let us now derive Result 2 for the pitchfork bifurcation. We expect that x^'^^ and 
C| are of order e and C„ = 0(6") from the stationary solutions. Thus we introduce 
the following multiple scales Ansatz: 

2 

i? = e^eJi?(i)(r,s) + 0(e^), (55a) 



i=o 






e^e^r(^)(r,s) + 0(e^ 



j=0 



CI 



(0) 



^eVi^-)(r,.) + 0(e^ 
Lj=o 



eV. 



(556) 



(55c) 
(55 (i) 



All the unknown functions in these equations should be bounded in the limit r — t- oo. 
Note that Lp^ = Lpf = 5jo. The Ansatz (155 cp follows the form of the equilibrium 
solution and the slowly varying time scale s is chosen to keep the functions i?'-^^ and 
r*^^) bounded as r — )• oo (see below). Note that (\21bf\ implies 

R 

1 3R — 777-7; — : 3i? 



X 



(0) 



c\ 



(0) 



a;(°) sin x 



X 



(0)3 



^(0)3_^^(0). 



X 



(0)3 



smx 



.(0)2 



(4Cr^ + i?^ 



(56 a) 
(566) 



These equations are needed to calculate 7 and 7 sin x up to terms of order e^. Inserting 
( 15^ and (I55ap -( 156B in (122 ap and l\22b\i and equating like powers of e, we obtain a 
hierarchy of equations. The leading order equations are: 

9,i?(°) = -i?(°\ 9,r(°) = 0, (57a) 

5.^i°^= 2(^11 -^f), ^r^l. (576) 

The solutions are i?(°) = R^°\0)e-^, cp^n^ = 1 + 0{e-^^) (for n > 2) and r(°) = r(°)(s). 
For large enough values of the fast time r, these solutions become 

i?(o) = o, r(°) = rW(s), </^i°) = i, (58) 

except for exponentially decreasing terms that we will ignore in what follows. We now 
find an equation for r*^°^(s). Equations (I56ap - fl566p become 

2 



j=0 



a:(°) sin x 



-(i?«+ei?(2)) 



X 



(0)3 



^e^r^")- 



X 



(0)3 



smx 



0. 



Then the equations for R^^^ and T^^^ with j = 1, 2 are 



2r(o)3 



(59 a) 
(596) 

(60a) 
(606) 
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Again R^^'^ = R^'^'^ = up to terms that decrease exponentially in the r scale. The 
solution r^^) of (16061) is bounded as r — )■ oo only if the right hand side of the equation 
vanishes, which produces the sought amplitude equation: 



rfrw 



20,rw - - r(o)=^ 



ds ^ 3 

The solution of (EH) for nonzero initial conditions is 



r(°)(s) = sign[r 



(0) 



\ 



3 02 



-46*2 5 



3 02 

r(o){o)2 



(61) 



(62) 



As the slow time s — ?■ oo, this function tends to if ^2 = —1 (d > 1) and to ±a/3 if 
02 = I {9 < 1). In the original variables, (|62|) becomes (131 a|) . These solutions tend 
to the stationary values C| ~ ±y^3(l — 0) corresponding to the stable equilibrium 
solutions if ^ < 1. For r(°)(0) = 0, it is T^^\s) = for all s. 



5. 3. Derivation of Result 3 

As ^ — ;■ 0, the reduced equations (122 gp and ( 122 6p for the correlations Cii and for the 
oscillation amplitude R can be simplified. Note that for x'-^-' 7^ the argument in (|20|) 
may become arbitrarily large, and therefore we have 



7 



[signxW](l-2e-^l^'°'l/^). 



If |C| I < i? we have that x'^^' changes sign due to the oscillatory term and 



7 ~ [sign C[ 



Xo 



arcsm 




.(0) 



(0), 



viR-\cr\)+vi\cr\-R) 



(63) 

(64a) 
(646) 



where ri{x) = 1 for x > and otherwise. We also find (see Appendix C ) 

2 cos Xo 2 



■2a;(o)' 
7 1 —^ ) smx 

7 I —;r- ) sm X 



TC 



IT 



2e-^l^'' 1/^/1 ( 




ioi \Cr\<R, (65a) 



(656) 



Using Eq. (06]) and ([63]) , we get 7 ~ ^ = signCf^ = ±1 for \c[°^\ > R, which is ( 132^ . 
The solution of the latter equation can be written as Cn = 5*" + A„(r), where A„ 
also satisfies the equations ( 132 ap but with boundary condition Aq = 0. This system of 
equations can be solved by using a generating function whose moments are the A^. We 



show in Appendix C that the result is 



CW(r) = S-- e-'^ 5^[5^' - Cf (0)][/|„_,|(2r) - /„+,(2r)] . 



(66) 
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where In{x) is the modified Bessel function of the first kind with integer index n. The 
correlations in ([MD decay algebraically dn\r) ~ 5" - n/(47rr3)i/2 Y,%ij[S^ - Cf\<S)\ 
as r — )■ oo (long time tails) [17]. In fact, this is the particularization of Glauber's 
solution for the Ising model p] to the zero temperature limit. 

According to (122^1) and (165 611 . R satisfies (132^1) . whose right hand side is 
exponentially small. Thus R{t) decays to zero only on an exponentially long time 
scale. Therefore C| approaches 5 = ±1 whereas x*^*^-* oscillates on the time scale r 
with almost constant amplitude -R(O) < 1 about it. It takes an exponentially long 
time for the oscillation amplitude to vanish. For equilibrium-like initial correlations 
C^°)(0) = r", with < r < 1, we show in |Appendix"C] that Ci°^(r) satisfies ^ 



For R > icf^l, §Ml and §M) yield §3^ and §M)- R decays until R ~ jcf^ 
Recapitulating, Cn and R are given by ( 132 aP and ( 132 6|) . respectively, if _R < |C| | and 
by ( 133 al) and ( 13361) if i? > |C} |. In any case, the correlations Ci tend algebraically to 
5" = (±1)" as Ci°^ - ^" oc r-3/2 (^ _^ oo)_ According to (l2T6l) . the oscillator vibrates 
with its natural frequency about an equilibrium position that approaches slowly a; = ±1, 
with an amplitude i? = 0(1), which vanishes only after an exponentially long time. 

6. Conclusions 

We have analysed the dynamics of a fast oscillator coupled to a one-dimensional chain 
of N Ising spins ai, i = 1, ... ,A^ in contact with a thermal bath at temperature 6. 
The interaction energy between the oscillator and the spins is proportional to the 
oscillator position and to ^^ctjCTj+i. In the limit of infinitely many spins, there is a 
second order phase transition at the critical temperature 6 = 1 and the fiuctuations of 
the oscillator position and the spin correlations are negligible. Provided the oscillator 
natural period is much smaller than the relaxation time of the spins, the oscillator 
position performs modulated oscillations about its slowly varying equilibrium which is 
related to the correlation between nearest neighbor spins. The spin correlations decay 
to their equilibrium values over the slow spin relaxation time scale r. Using a multiple 
scale analysis that exploits the large separation between characteristic times, we derive 
modulation equations for the average oscillation envelope and phase of the oscillator 
and for the spin correlations. 

The modulation equations can be analyzed in two limits: near the critical 
temperature 6 = 1 and in the limit of low temperatures. At the critical temperature, the 
"trivial" solution of the modulation equations (zero spin correlation and zero oscillator 
position) undergoes a supercritical pitchfork bifurcation: the trivial solution is stable 
for 6 > 1 and two non-trivial stable stationary solutions bifurcate from it for 6 < 1. We 
have constructed these solutions using multiple scales analysis and shown the typical 
slowing down near the critical temperature. Below the critical temperature, a study of 
the modulation equations show that most initial conditions evolve towards one of the 
stable stationary solutions revealed by the bifurcation analysis. The spin correlation 
Ci = (aj(Tj+i) and the average oscillator position corresponding to these stationary 
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solutions both approach ±1 as — )■ 0. However, the instantaneous spin correlation for 
a fixed small temperature tends to a stable stationary value algebraically as r — )■ oo. 
The oscillator position carries out modulated oscillations about the instantaneous spin 
correlation with an amplitude that decays to zero over an exponentially long time. 

In the thermodynamic limit A^ ^ 1, perhaps the most striking result is the 
dynamical stabilization of the unstable trivial equilibrium state Xeq = C*i,eq = at any 
temperature 6 < 1 for random initial conditions of the Ising system. In the limit of a 
very fast oscillator we are considering, the analysis of the averaged modulation equations 
shows that the oscillator position averaged over one natural period, x(t), tends to 
approach its instantaneous rest position Ci{t) (which evolves much more slowly towards 
equilibrium). But, for random initial conditions of the spin system Ci{t) = 0, and the 
envelope of the oscillator position tends to zero with relaxation time 6/6. We have 
verified our theoretical predictions by direct numerical simulation of the oscillator-spin 
system and by numerically solving the macroscopic equations which ignore fiuctuations 
and are valid in the thermodynamic limit. In general, the agreement between simulation 
and theory is excellent. For large enough size A^ and random initial conditions, direct 
simulations show that the fiuctuations are so small that the system is not able to depart 
from the unstable state, in every single trajectory. For smaller systems, fiuctuations are 
larger and drive both the oscillator and the Ising system to one of the two stable states. 

These results are somewhat surprising given that there is a Lyapunov functional 
for the oscillator-spin system guaranteeing that any initial condition evolves on the 
time scale r towards a stable stationary state corresponding to a canonical equilibrium 
probability distribution. However, the corresponding H-theorem is proved for finitely 
many spins and therefore it does not necessarily applies to fiuctuation-free results 
obtained after the limit A^ — )■ oo has been taken. Stabilization of the unstable trivial 
stationary state for random initial conditions is reminiscent of weak ergodicity breaking 
in glasses p3| [T6] . It is different from noise- induced stabilization and selection of 
unstable states as studied by Freidlin [18] and Muratov et al |T9]. 
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Appendix A. H-theorem 

Let us define the functional 

i/w^/d./dp5:n..P.<..«)i.(|^^) . (A.1) 
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where V{x,p,cr,t) and Vo{x,p,(T,t) are two different solutions of the master equation 
([6]). We have H{t) > because x Inx > x — 1 for all x = V/Vq > and the probability 
densities are normalized to 1. By using ([6]) and integrating by parts, the time derivative 
oiH, 

(V 



Hit) = dx dpY^ 

can be shown to be 

H{t) = S I dx 



dfVln 



V_ 



[v. ) '^^^ 



(A.2) 



'^pJ2J2v'' 



V{x,p,cr,t) 



Vo{x,p,cr,t) 
X [W{cr' -)■ cr\x,p)V{x,p,cr',t) - W{cr -> a'\x,p)V{x,p,cr,t)] 

— [W{a' — > a\x,p)Vo{x,p, cr' ,t) — W{a — > a'\x,p)VQ{x,p,cr,t)] 
V{x,p,(T,t) 



X 



Vo{x,p,(T,t) 

Here we have used the notation 
W{(r -> (r'\x,p) = 



Wi{cr\x,p) if a' = RiCr, i 
otherwise, 



1 N 



(A.3) 



(A.4) 



where RiCr is the configuration obtained from cr by flipping the i-th spin. By 
interchanging cr and <t' in the second term on the rhs of flA.Sp . and defining 

V{x,p,a,t) 



f{x,p,cr,t) 

Vo{x,p,a,t) 

equation (\A.3\i can be rewritten as 

H{t) = 6 dx dp^Y^ W{cr' -^ a\x, p) {P(x, p, cr', t) 

X [ln/(x,p, cr,t) -ln/(x,p,cr',t)] 

- 'Po(a;,p,cr',t)[/(x,p, cr,t) - /(x,p, cr', t)]} . 

Let / = f{x,p,(T,t) and /' = f{x,p,(T',t). Then we can write (JA.6P as 



(A.5) 



(A.6) 



/ + / 



-6 



X 



(A.7) 



f)Kf)-(f)-l-- 

Since W fV^ > and x In x — x+1 > for x > 0, -ff (t) decreases with time. Moreover the 
expression in square brackets in flA.7p vanishes only if / = /'. Thus H{t) > decreases 
monotonically as time increases and is bounded from below. Hence H{t) tends to a 
limit as t — )• oo, such that Imit^ooHit) = 0. According to (lA.7p . this is possible only 
if / = /' for all pairs of connected states having W^cr' -^ cr|x,p) 7^ 0. For fixed x,p, 
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all the spin configurations a are connected through a chain of transitions with non-zero 
probability, and therefore 

lim r';'"-"-^' =/fep). (A.8) 

independently of the spin configuration a. We have shown that in the long time limit, 
the ratio of any two solutions of the master equation is a function of x and p. Since 
the equilibrium distribution Veq{x,p,a) is also a solution of the master equation, any 
time-dependent solution satisfies 

V{x,p,(T,oo) =Vcciix,p,(T)f{x,p). (A.9) 

The distribution V{x,p, a, oo) must be a stationary solution of the master equation IQ. 
By inserting (I A. 9 1) in (E]) and using detailed balance, we find that f{x,p) is constant, 
independent of x and p. Clearly f{x,p) = 1 due to the normalization condition for 
V{x,p,cr, oo). Then 

V{x, p, (T, oo) =Vcciix,p, cr) and hm ' — '-— = I. (A. 10) 

t^ooVo{x,p,(T,t) 

Equation ( ]A.10|) establishes that all the solutions of the master equation corresponding 
to N spins tend to the canonical equilibrium distribution in the long-time limit. As 
N ^ oo, there is a second order transition at the critical temperature 6 = 1. For 6 > 1, 
there is only one phase given by the limit of the canonical equilibrium distribution as 
N ^ oo with order parameter x^q = Ci^q = 0. For 9 < 1, there are three phases 
which appear as limits of the canonical distribution with order parameters Xeq = 
(unstable phase) and ±|xoq| 7^ (stable phases). These phases coalesce at 6* = 1 
in a pitchfork bifurcation, as described in [10]. Each sign attracts different initial 
conditions. The multiplicity of (macroscopic) equilibrium solutions corresponding to 
extrema of the equilibrium distribution is a direct consequence of the nonlinearity of 
the macroscopic equation, which is compatible with the linearity of the master equation 
and the multiplicity of its stationary solutions as iV — )■ oo [20] . 



Appendix B. Pitchfork bifurcation in the macroscopic equations 

Let us assume that 6 = 0(1) and seek a solution of the macroscopic equations (TT8l)-( l20|) 
near the critical temperature 6 = 1 by means of the multiple scales Ansatz: 

2 

x = eJ2e'^'^'\t,0 + 0{e'), (B.l) 

3=0 

2 

C; = e^e^r(^)(t,C) + 0(e^), (B.2) 



Cn — ( C*l 



2 



.4\ 



(B.3) 



^eVi^-)(t, + 0(6^ 

.3=0 

9 = 1- 92e\ (B.4) 

C = eh. (B.5) 
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As in section El e measures the amplitude of the bifurcating solutions and all unknowns 
in f lRT]) - flRSJ) should be bounded as t ^ oo. Inserting flRTl) - flRSl) in ([l8])-([20]) and 
equating like powers of e, we get the following hierarchy of equations 

(^ + l) e^ - r(°) = 0, (B.6) 

(^^ + 26\ rW - 25e(°) = 0, (B.7) 

(^ + l) ^(^^ - r(^) = 0, (B.8) 

(^ + 25\ r(i) - 2(5^(1) = 0, (B.9) 

f^ + 25") (rwv5°^) = 25^(°)rw, (B.12) 

and so on. 

The solution of flRGjl - flRT]) is 

^(0) = K{C) + Re [A{C)e^'^^^^'] , fi = VT^^, (B.13) 

rW = ir(C) - 26 Re [A(C)e(*^-'^)*(«n - 6)] , (B.14) 

for < 5 < 1 and similar formulas for the case S > 1. In both cases, we can ignore 
terms that decrease exponentially rapidly in the fast scale t and set 

^(0) = K{C) + EDT, r(°) = K{C) + EDT, (B.15) 

where EDT stand for exponentially decreasing terms in the fast time scale. To find 
K{C), we insert flRTSl) in fiRTOl) - f lRT2l) . ignore EDT and obtain: 

(^ + 1) e^^^ - r(^) = 0, (B.16) 

(I + 25] r(^) - 2<5e(^) = 25K U + r(°)Vf - IkA - ^, (B.17) 
f^^ + 2(5') (r(o)Vr) = '2SK\ (B.18) 

The solutions of (IRIHI) and ( iBlSJ) are ^(2) = r(2)+ EDT and 1(0)^2°^ = K^+ EDT, 
respectively. Then (]B.17p becomes 

9r'^' 2«f«,-M-^, (B.19) 



dt \ 3 J dC 

plus EDT which have been ignored. The solution F*^^) is bounded as t — )■ 00 provided 

f^2.A-(<,.-lA-). (B.20) 
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This amplitude equation is the same as (1611) and therefore it has the same solution 

K{eH) = [signir(O)] 



+ 0(|1-^|). (B.21) 



^ 1 + e-^i^-m (Ml_^ _ i) 



In the original variables, we recover x = Ci with Ci given by fl31ap . i.e., Result 2 is 
valid both for the macroscopic equations and for their averaged version in Result 1. If 
we keep EDT in the leading order approximation, we obtain the following composite 
expansion instead of fl31ap - fl31cp : 

x{t) ~ K{eh) + a(0) e~^^ cos(fit + ^), (B.22) 

Ci{t) ~ K{eH) + 2(5 a(0) e-^'[6 cos(fit + if) + n sin(fit + ^)], (B.23) 

K{0) = r + 26vo, Ci{0) = r, x(0) = xq, ^(0) = vo, (B.24) 



a(0) = y (.0 - r - 26voy + ^l - 26^H S^ix, - r)]^ ^ ^^_^^^ 

Xn — r — 26Vr) ,^ 

Note that K{0) — )■ r and that a(0) and y? become -R(O) and 7r/2 — given by ( 1261) as 
5 — )• 0. The approximation (1B.22P -( 1B.26P is better than f l31ap - fl31d) for relatively large 



values of 6 but both approximations become indistinguishable as 5 — >■ 0. 

Appendix C. Averages at low temperature 

Let us assume that |CJ | < R and that CJ < for example. Let us consider the 
limit e -^ 0+. Then a;(°) = C[°^ + Rsinx > and therefore 7(2x(°V^) ~ 1 if 
Xo<X — ^ + 0<^~Xo (sinxo = IC*} \/R with < xo < 7r/2 as in 164 ^jP and 
x^^^ < 0, 7(2x*^°V^) ~ ~1 fo^ ~^ < X < Xo or TT — Xo < X < TT. Then 



and the other possibilities in (164 gp are obtained using similar calculations. 
For < -C{°^ < R, we find 



7 ( —g- ) ^^^ ^ ~ 27r I / ^^^ "^ ~ / ^^^ "^ ~ / ^^"^ 

1 r-^" . 2COSX0 .^^. 

= -/ smxdx= , (C.2) 

which is ( 165 gp (the case < C| < -R gives the same formula). Eq. ( 165 6p follows from 
using ln(l + x) ~ x in ( 152|) and the integral formulas for the modified Bessel function. 
For CJ > -R, the equations for r„ = Cn — 1 are 

-^ = r„_i + r„+i - 2r„, Fq = o. (c.3) 

dr 
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If we define r_„ = — r„, the generating functional 

oo oo 

r(g,r) = Y. e^"'r„(r) = 2z^r„(r) sinM (C.4) 



n=l 



obeys the equation 
dV 



= _2(1 - cosg)r, r(g,0) = 2^ J][C(°)(0) - 1] sin(ng). 



dr 



n=l 



The solution of (IC.SP is T(q, r) = T{q, 0) e "^^^ cosg)r^ from which we get 
CW(r) = 1 - - rT{q,0)e-^^'-'°"'^^sm{nq)dq. 



(C.5) 



(C.6) 
■" Jo 

Inserting the initial condition (JC.SP in ( ]C.6P and using the integral representation 

formula for the modified Bessel functions [17], we obtain ( l66ll with 5 = 1: 

oo 

Cf (t) = 1 - e-^^5^[l - Cf (0)] [/|._,|(2r) - /.+,(2r)]. (C.7) 

i=i 

On the other hand, for an equilibrium-like initial condition, Cn (0) = r", < r < 1, 
the initial condition in (IC.5I) becomes 



r(g,0)=2« 



r smg 



smg 



1 + r^ — 2rcosg 2(1 — cos g) 
Inserting ( 1C.8|) in ( ]C.6|) . we obtain 



1 f^ sing sin (ng) 



Qn{r) : , ^ 

TT Jo 1 ~ COS g 

Differentiating (IC.lOll . we find 



^_^ 2 r rsingsin(ng) ^_2(i-cosa)T 
TT Jo 1 + r^ — 2r cos g 



+ Qn{r) 



-2(1— cos q)T 



dq. 



(C.8) 

(C.9) 
(CIO) 



dr 



TX 



sing sin (ng)e' 



2(1— cos(jr)r 



dq 



,-2t /"vr 



vr 

-2t| 



/ {cos[(n - l)g] - cos[(n + \)q\}e^^^'''\ 
Jo 



whereas 



Qn(0) 



e-'"[/„_i(2r)-/„+i(2r)], 
1 [^ sing sin (ng) 



TT 



1 — COS g 



(ig 



The solution of (]C.11|) with initial condition f lC.12p is 

Qn{r)= [ e-2*[/„_i(2t)-/„,+i(2t)]rft-l. 
Jo 

Equations (109]) and ( lal3|) yield (|38|). 



(C.ll) 
(C.12) 

(C.13) 
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